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Introduction 



Recent years have brought about exciting new developments in computerized 
tomography. In particular, a novel, very promising approach to the creation 
of diagnostic techniques consists in combining different imaging modalities, 
in order to take advantage of their individual strengths. Perhaps, the most 
successful example of such a combination is the Thermoacoustic Tomog- 
raphy (TAT) (also called photoacoustic tomography and optoacoustic to- 
mography and abbreviated as TCT, PAT, or OAT) [1-8]. 

Major progress has been made recently in developing the mathematical 
foundations of TAT, including proving uniqueness of reconstruction, obtain- 
ing range descriptions for the relevant operators, deriving inversion formulas 
and algorithms, understanding solutions of incomplete data problems, sta- 
bility of solutions, etc. One can find a survey of these results and exten- 
sive bibliography in [9]. In the present article we concentrate on the recent 
advances in the inversion formulas and algorithms for TAT. Mathematical 
problems of the same type arise also in sonar, radar, and geophysics applica- 
tions (e.g., [10-12]). Discussion of some mathematical problems concerning 
TAT can be also found in the chapters written by D. Finch and Rakesh and 
by S. Patch. 

While this text addresses the mathematics of TAT only, one can find 
extensive discussion of physics, engineering, and biological issues related to 
TAT in the recent surveys [4,5,8], textbook [7], as well as in other chapters 
of this volume. 
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1 Thermoacoustic tomography 



We give first a brief description of TAT. The data acquisition starts with a 
short electromagnetic (EM) pulse being sent through the biological object 
under investigation (e.g., woman's breast in mammography) . A fraction of 



Transducer 




MW pulse 



Figure 1: The TAT procedure. 

EM energy is absorbed at each location x inside the object, thus triggering 
thermoelastic expansion of the tissue and emergence of a pressure wave p(x, t) 
(an ultrasound signal) that, in turn, is measured by transducers placed along 
some observation surface S surrounding (completely or partially) the object. 
The initial pressure po(x) = p(x,0) is determined by the intensity of the 
EM pulse (that assumed to be known) and by the local properties of the 
tissue. It is known (e.g., [1,4,5,8, 13]) that in the radiofrequency and visible 
light ranges absorption of the EM energy by cancerous cells is several times 
stronger than by the healthy ones. Thus, knowledge of the initial pressure 
Po(x) would provide an efficient tool for early detection of cancer. Frequently, 
the ultrasound contrast is sufficiently small to justify the use of the constant 
sound speed approximation. Most work on TAT up to date is based on 
this assumption. However, such an approximation is not always appropriate; 
some of the results described below, as well as in [9, 14, 15] aim towards the 
general case of a variable speed of sound. 

Once the data p(x,t) has been measured on S x IR + , one can attempt 

1 It has been argued that the radiofrequency and visible light ranges are most appro- 
priate in TAT [8]. For the purpose of this text, no distinction is made between these 
cases. 
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to recover from p(x,t) the initial value po(x) of the pressure inside S (the 
thermoacoustic image). 



2 Mathematical model of TAT 



Let us for notational convenience denote po(x) (the image to be recon- 
structed) by f(x). In this section, we present a mathematical description 
of the relation between the functions f(x) and p(x,t). We assume that the 
function f(x) is compactly supported in R n (we allow the dimension to be 
arbitrary, albeit the most interesting cases for TAT are n = 3 and n = 2). 
At each point y of an observation surface S one places a point detector 2 
that measures the value of the pressure p(y,t) at any moment t > 0. It is 
usually assumed that the surface S is closed (rather than, say, cylinder or a 
plane 3 ). It is also assumed that the object (and thus the support of f(x)) is 
completely surrounded by S. The latter assumption is crucial for the validity 
of most inversion formulas; however in some cases we will be able to abandon 
this requirement. 

The mathematical model described below relies upon some physical as- 
sumptions on the measurement process, which we will not describe here. The 
reader can find such a discussion in [8]. 

We assume that the ultrasound speed v s (x) is known, e.g., through trans- 
mission ultrasound measurements [15]. Then, the pressure wave p(x,t) sat- 
isfies the following set of equations [13,23,24]: 



Now one needs to recover the initial value f(x) at t — of the solution p(x, t) 
from the measured data g(y,t) := p(y,t),y e S, t > 0. Incorporating this 

2 Planar and linear detectors have been considered as well, see [16, 17] and further 
references in [9]. 

Reconstruction formulas for the planar and cylindrical cases are well known, see e.g. 




n 



(1) 



[18-22]. 
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data, one rewrites (1) as 



'ptt = v s 2 {x)A x p, t>0, x G W 
^ p(x,0) = f(x), 
* Pt{x,0)=0 

^p{y,t)=g{y,t), yeSxR + 





Figure 2: An illustration to (2). 

In other words, we would like to recover the initial value f(x) in (2) from 
the knowledge of the lateral data g(y,t) (see Figure 2). At a first glance, it 
seems that the data is insufficient for the reconstruction, i.e. for recovering 
the solution of the wave equation in a cylinder from the lateral values alone. 
However, this impression is incorrect, since there is additional information 
that the solution holds in the whole space, not just inside the cylinder Sxl + . 
To put it differently, if one solves not only the internal, but also the external 
problem for the wave equation with the data g on the cylinder S x IR + , 
then the solutions must have matching normal derivatives on S x M + . In 
most cases, this additional information provides uniqueness of recovery of 
f(x) (see below, as well as [9,14,25-29], and references therein). It is also 
sometimes useful to notice that p can be extended as an even function of 
time and thus satisfies the wave equation for all values of t. Similarly, data 
g can be extended to an even function. This, in particular enables one to 
apply Fourier transform in time. 

An additional structure arises in this problem, if one assumes that the 
object under investigation is nearly homogeneous with respect to ultrasound: 
v s (x) = 1- In this constant speed case, there is an alternative way to describe 
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the relation between the data g(y, t), (y, t) G S x R + and the unknown image 
f(x),x G M 3 . The known Poisson-Kirchhoff formulas [30, Ch. VI, Section 
13.2, Formula (15)] for the solution of (1) with v s — 1 give 

p(x,t) = ^(t(Rf)(x,t)), (3) 

where 

(Rf)( x ,r) = j- J f(x + ry)dA(y) (4) 
li/l=i 

is the spherical mean operator applied to the function /(x), and dA is the 
surface area element on the unit sphere in IR 3 . Thus, the function g(y,t) 
for y G 5 and all t > essentially carries the same information as the 
spherical mean Rf(y : t) at all points (y,t) G S x M + (see, e.g., [27]). One 
can, therefore, study the spherical mean operator R : f — > i?/ and, in 
particular, its restriction i?s to the points y G 5 of the observation surface: 

R s f(x,t)= J f(x + ty)dA(y), x e S, t>0. (5) 
l»l=i 

This explains why in many studies on thermoacoustic tomography, the spher- 
ical mean operator has been used as the model. One needs to notice, though, 
that in the case of a non-constant sound speed, the spherical mean interpre- 
tation (as well as any integral geometry approximation) is no longer valid, 
while the wave equation model still is. 



3 Uniqueness of reconstruction 

Uniqueness of reconstruction of a compactly supported (or sufficiently fast 
decaying) function f(x) from the data g collected from a closed surface S is 
well known in the case of a constant sound speed (i.e., when the interpretation 
in terms of spherical mean operators is possible). One can find discussion of 
such results in [9, 14, 25, 27-29, 31-34]. 

In the case of a variable sound speed, it is shown in [31, Theorem 4] that 
uniqueness of reconstruction also holds for a smoothly varying (positive) 
sound speed, if the function f(x) is supported inside the observation surface 
S. The proof uses the famous unique continuation theorem by D. Tataru 
[35]. 
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We present now a recent simple uniqueness theorem that also allows a 
non-constant sound speed v s (x) and does not require the function to be 
supported inside S. In order to do so, we need to formulate first some 
assumptions on v s (x) and the function f(x) to be reconstructed. 

1. Support of f(x) E #f oc (R n ), s > 1/2 is compact. 

2. The sound speed is smooth (a condition that can be reduced), strictly- 
positive v s (x) > v o > and such that v s (x) — 1 has compact support, 
i.e. v s (x) = 1 for large x. 

3. Consider the Hamiltonian system in with the Hamiltonian H = 



The solutions of this system are called bicharacteristics and their pro- 
jections into are rays. 

We will assume that the non-trapping condition holds, i.e. that all 
rays (with £ ^ 0) tend to infinity when t — > oo. 

Theorem 1. [14] Under the assumptions formulated above, compactly sup- 
ported function f(x) is uniquely determined by the data g. (No assumption 
of f being supported inside S is imposed.) 

Uniqueness fails, however, if / does not decay sufficiently fast (see [25], 
where it is shown for the constant speed in which spaces L p (R d ) of functions 
f(x) closed surfaces remain uniqueness sets). 

4 Reconstruction in the case of constant sound 
speed: formulas, algorithms, and examples. 

We consider here the case of a constant sound speed: v s (x) = 1. One can 
work then either with the wave equation, or with the spherical mean operator 
model. 
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4.1 Inversion formulas and procedures 

Consider the case of the observation surface S being a sphere. The first 
inversion procedures for this situation were obtained in [36] in 2D and in [37] 
in 3D by harmonic decomposition of the measured data g and of the function 
/, and then by equating coefficients of the corresponding Fourier series (see 
also [9] for a brief description of this procedure). The two resulting series 
solutions are not quite analogous. Indeed, in [36] one had to divide the 
Hankel transform of the data by the Bessel functions that have infinitely 
many zeros, which would create instabilities during implementation. The 
3D solution in [37] is free of this difficulty and can also be adopted for 2D. 
We will see a different type of series solutions later on in this section. 

4.1.1 Approximate inversion formulas 

The standard way of inverting Radon transform in tomographic applications 
is by using filtered backprojection type formulas [20,38-41]. It combines a 
linear filtration of projections (either in Fourier domain, or by a convolution 
with a certain kernel) followed (or preceded) by a backprojection. In the case 
of the set of spheres centered on a closed surface (e.g., sphere) S, one expects 
such a formula to involve a filtration with respect to the radial variable and 
an integration over the set of spheres passing through the point x of interest. 
Albeit for quite a long time no such formula had been discovered, this did 
not prevent practitioners from reconstructions. The reason was that good 
approximate inversion formulas (parametrices) could be developed, followed 
by an optional iterative improvement of the reconstruction [6, 13,21,22,42- 
44]. 

Perhaps the most advanced approach of this kind was adopted by Popov 
and Sushko [42,43]. These authors have developed a set of "straightening" 
formulas that allow one to reconstruct from the spherical means an approxi- 
mation to the regular Radon projections. The main idea is that for each (hy- 
per)plane passing through the support of the function to be reconstructed, 
one builds a family of spheres with centers at the detectors' locations and 
tangential to that plane. One such sphere is chosen for each point of the plane 
contained within the support. The integrals over these spheres are known, 
as they form a subset of projections g. An approximation to the integral of 
the function over the plane is then computed by integrating over these pro- 
jections a functional (local in odd and non-local in even dimensions). When 
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all the plane integrals are computed, the function is reconstructed by apply- 
ing inversion formulas for the regular Radon transform. This procedure is 
not exact; however, as shown in [42], such an algorithm yields a parametrix. 
Namely, the difference between such an approximation and the original func- 
tion / is described by a pseudodifferential operator of order —1 applied to 
/. In other words, reconstruction is accurate up to a smoothing operator. 
This result holds even if the measuring surface is not closed (but satisfies a 
"visibility" condition), which is important for applications in the problems 
with incomplete data. 



4.1.2 Exact filtered backprojection formulas in 3D 



The first set of exact inversion formulas of the filtered backprojection type for 
the spherical surface S was discovered in [29]. These formulas were obtained 
only in odd dimensions (and then extended to even dimensions in [45]). Var- 
ious versions of such formulas (different in terms of the order in which the 
filtration and backprojection steps are performed) were developed. 

To describe these formulas, let us assume that B is the unit ball, S = dB 
is the unit sphere in M 3 , and a function f(x) is supported inside S. The 
values of its spherical integrals g(z,r) with the centers on S are assumed to 
be known: 



g(z,r) = J f(z + rs)r 2 dA(s) = Anr 2 R s f(z,r), z e S. 



(7) 



Some of the 3D inversion formulas of [29] are: 

g(z, \z-y\) 



f(y) 



"8^ 



dA(z), 



(8) 



f(y) 



8tt 2 



1<P_ 
tdt 2 



g(z,t) 



dA(z). 



t=\z-y\ 



(9) 



A different set of explicit inversion formulas, which work in arbitrary 
dimensions, was found in [46]. In 3D case the general expression derived 
in [46] simplifies to 



1 d g(z,t) 
tdt t 



dA(z), 



(10) 



t=\z-y\ 
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where n(z) is the vector of exterior normal to S. (We eliminated in this ex- 
pression the minus sign erroneously present in the original formula.) Equa- 
tion (10) is equivalent to one of the 3D formulas derived earlier in [47]. 

Similarly to the case of the standard "flat" Radon transform, all these 
3D inversion formulas are local, i.e. in order to reconstruct a value of the 
function at a certain point, one needs to know only values of all the integrals 
over the spheres passing through an infinitesimally small neighborhood of 
that point. 

It is worth noting that although formulas (9) and (10) yield identical re- 
sults when applied to functions that belongs to the range of the spherical 
mean Radon transform, they are in general not equivalent, i.e. lead to differ- 
ent reconstructions when the data is outside of the range (for instance, due 
to errors). Another important fact about these reconstruction techniques is 
that, unfortunately, they do not yield correct reconstruction within the re- 
gion surrounded by the detectors if the source is not contained within this 
region. Both these statements can be easily proven by the following example. 
Let us assume that the source function f(x) is constant (equal to 1) within 
the ball -8(0, 3) of radius 3 centered at the origin. In order to reconstruct the 
function within the unit ball, both formulas (10) and (9) use only integrals 
over spheres with the radius less or equal to 2, and centered at the points of 
the unit sphere. Obviously, all these spheres lie within the B(0, 3), and thus 
the projections g(z, t) are equal to the areas of the corresponding integration 
spheres, i.e. to 4nt 2 . By substituting this expression into (9), we obtain 



Function fi(y) defined by the above formula is harmonic in the interior of B, 
since the integrand is the free space Green's function of the Laplace equation. 
Due to the symmetry of the geometry, fi(y) is radially symmetric, i.e. it 
depends only on \y\. Therefore fi(y) = const for all y e B \ S. Let us 
compute /i(0): 




s 




s 



Thus, h{y) 



4 for all y e B \ S. 
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A similar computation with the use of (10) yields 



If 1 

/ 2 (y) = — div / n(z)- -dA(z) 

Zn J \z - y\ 



d 1 , x 47T 

• ,— m = 7; — = 2, 

2tt 7 dn(z) |^ — 2/| w 2tt 

s 

where we used the 3D Gauss formula. Both results f\ and fi are incorrect 
(not equal to 1). Besides, they are different, which proves that formulas (9) 
and (10) are not equivalent. 

One of the important benefits of having exact inversion formulas is that 
often a rather straightforward discretization of such a formula yields an ef- 
ficient and stable reconstruction algorithm. Such algorithms were developed 
in [48] using equations (8) and (9), and in [46] utilizing formula (10). 

In the simplest case, when the image is reconstructed on a grid of size 
m x m x m from 0{m 2 ) projections, each of which contains values for 0(m) 
integration spheres, all these algorithms have complexity of 0(m 5 ) opera- 
tions. In practical terms, for m of order of a hundred, the reconstruction 
time is measured in hours. An example of the reconstruction in 3-D using 
a method based on formula (10) is shown in Fig. 3. Reconstructions us- 
ing formulas (8) or (9) are quite similar in terms of stability, accuracy, and 
computation time. 





Figure 3: A mathematical phantom in 3D (left) and its reconstruction using 
inversion formula (10) 
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4.1.3 Exact filtered backprojection formulas in 2D 

Exact inversion formulas were obtained for even dimensions in [45]. Denot- 
ing by g, as before, the spherical integrals (rather than averages) of /, the 
formulas in 2D look as follows: 



2R 



Att 2 R 




A / / g(z,t)\og\t 2 - \y - z\ 2 \ dt dl(z), 



(11) 



S 



or 



2R 



m = 



4tt 2 R 




d (. d g(z,t) 



T7T t 



S 



dt V dt t 



log |t 2 — \y-z\ 2 \ dtdl(z), (12) 



where B is a disk of radius R centered at the origin, and S = dB is its 
boundary. 

Another 2D inversion formula [46] takes the following form (again, cor- 
rected for a sign): 



f(v) = ~g~ div J n ( z )K z , \V - z\)dl(z), 



(13) 



where 



h(z,t) 



2R 



Yo(Xt) / J {\t')g{z,t')dt' 



2R 



J (Xt) I J Y (Xt')g(z,t')dt' 



XdX, 



(14) 



and Jo(t) and Y (t) are the Bessel and Neumann functions of order 0. By 
analyzing the large argument asymptotics of these functions one can see [46] 
that the filtration operator given by equation (14) is an analog of the Hilbert 
transform. 

This reconstruction procedure can be re- written in a form similar to (11) 
or (12). Indeed, by slightly modifying the original derivation of (13), (14), 
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one can obtain a formula that would reconstruct a smoothed version f(x, v) 
of f(x) defined by the formula 

f(x,v)=F- 1 (\Z\-Tf), 0<u<l, 

where J 7 , T~ x are correspondingly the 2D Fourier and inverse Fourier trans- 
forms. The restriction of f(x, is) to the interior of the disk B is recovered by 
the formula 

f(y, v) = ~^ div J n(z)h v (z, \y - z\)dl(z), (15) 
s 

where 

K{z,t) = J Y (Xt) j J J (\t')g{z,t')dt'\-J (\t) I J Y G (\t')g(z,t')dt' \ \~ v d\. 

R+ \0 / \0 / 

(16) 

For < is < 1, one can change the order of integration in (16) to obtain 

2R 

h u (z,t) = J g(z,t')K u (z,t,t')dt', (17) 
o 

K u (z,t,t') = j \ (\t)J (\t')\- u d\- J J (\t)Y (\t')\- u d\. (18) 

Using [49, formula 4.5, p. 211], the integral / Y (\t)J (\t')\- u d\ can be 
integrated exactly, yielding 



J Y (\t)J (\t')\- u d\ = 



^r(i-i/)j£^, t>t> 



-)1 — IS 



-ra-^^z^, t<t> 



The expression for the second integral in (18) is derived by interchanging t 
and t', which results in the formula 
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K u (z,t,t') = 



^{i-v y-'pyp?-" , t>t> 
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Finally, we substitute the above expression for K u (z,t,t') into (17) and take 
the limit v — > 0, to arrive at the following formulas 

f(y) = 7^ div J n(z)ho(z, \y - z\)dl(z), 
s 

2R 

h (z,t) = / g(z,t') 2 _ dt! 



or 



2R 



f(y) = ^div J n(z) J g(z,t') t/2 _l y _ zl2 dt' 
s Lo 



dl(z). (19) 



Similarly to the one appearing in (11) and (12), the filtration operator 
in (19) also involves kernel t ,i_ t i ■ If desired, it can be re-written in the form 
of a convolution, either by a change of variables t 2 — > t, or by noticing that 

2 l/t' l/t' 



t' 2 -t 2 t + t' t-t r 

This is important from the computational point of view, since it allows the 
reduction of the inner integral in (19) to the sum of two Hilbert transforms, 
computational algorithms for which are well known. 

All inversion formulas presented in this section require 0(m 3 ) operations 
to reconstruct an image on a grid of size mxm from 0{m) projections, each 
consisting of 0(m) values of circular integrals. This coincides with the oper- 
ation count required by a classical (non-accelerated) filtered backprojection 
algorithm in 2D. 

It is not yet known currently whether formula (19) is equivalent to (11) 
and (12). However, as shown in the previous section, this is not the case 
for the 3D versions of these formulas, and thus this seems unlikely in the 
two-dimensional well. 

Finally, similarly to the filtered backprojection formulas for the classical 
2D Radon transform, the inversion formulas (11), (12), and (19) are not 
local. In other words, in order to recover the value of f(x) for a fixed point 
x, all the values of g(z,t) have to be known. 
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4.2 Series solutions for arbitrary geometries 

Explicit inversion formulas for closed surfaces S different from spheres have 
not yet been found 4 , except the result of [14] described in the next Section. 
There is, however, a different approach [50] that theoretically works for any 
closed S and that is practically useful when the surface is the boundary of a 
region, in which the spectrum and eigenf unctions of the Dirichlet Laplacian 
are known (or could be effectively approximated numerically). 

Let (where \ k > 0) and Uk(x) be the eigenvalues and normalized eigen- 
functions of the Dirichlet Laplacian -A^ on the interior Q of the observation 
surface S: 



Au k {x) + X 2 k u k {x) =0, x en, flCR", (20) 
u k (x) =0, x E S = dil, 



KII2 = J \u k (x)\ 2 dx = 1. 



As before, we would like to reconstruct a compactly supported function f(x) 
from the known values of its spherical integrals g(z,r) (7). 

According to [50], if f(x) is represented as the sum of the Fourier series 

00 

f(x) = ^a k u k (x), (21) 

m=0 

the Fourier coefficients a k can be reconstructed as follows: 

f d 

a k = I(z,\ k )—u k (z)dA(z) (22) 
J an on 



where 



I(z,X k ) = J g(z,r)$ Xk (r)dr, 



and Q\ k (\x — z\) is a free-space rotationally invariant Green's function of the 
Helmholtz equation (20). 



4 Planar and cylindrical observation surfaces, for which such formulas are known [18-22], 
are not closed. 
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Formula (22) is obtained by substituting the Helmholtz representation 
for Uk(x) 



into the expression for the projections g(z,t). 

This eigenfunction expansion approach requires the knowledge of the 
spectrum and eigenfunctions of the Dirichlet Laplacian, which is available 
only for some simple domains. However, when this information is available, 
the method yields reliable, robust, and, in some cases, fast reconstruction. 
For example, as it was shown in [50], for the cubic observation surface S, one 
can compute reconstructions thousands times faster than by methods based 
on explicit inversion formulas of backprojection type discussed above. The 
operation count for such an algorithm is 0{rr? logm), as compared to 0(m 5 ) 
for the explicit inversion formulas. 

Another advantage of the series technique is its ability to "tune out" the 
signal coming from outside of S. In other words, unlike the explicit inversion 
formulas discussed in the previous sections, the present method enables one 
to reconstruct the values of f(x) for all x lying inside S even in the presence of 
the sources outside. We illustrate this property by the reconstruction shown 
in Fig. 4. (The dashed line in the left figure represents surface S, i.e., the 
location of the detectors.) 

5 Reconstruction in the case of variable sound 
speed. 

In this section we consider a more general case of the variable sound speed 
v s (x). Our analysis is valid under previously imposed conditions on this 
speed, namely, that v s (x) is sufficiently smooth, strictly positive, non-trapping, 
and v s (x) — 1 is compactly supported. 

Consider the Hilbert space H = L 2 (VL,v s ~ 2 (x)dx) ) i.e., the weighted L 2 
space with the weight v s ~ 2 (x). In this space, the naturally defined operator 



in f2 with zero Dirichlet conditions on S is self- adjoint, positive, and has 
discrete spectrum {\j.}(\k > 0) with eigenfunctions ipk( x ) G H. 




x G f2, 



(23) 



A = -v s 2 (x)A 
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Figure 4: The phantom shown on the left includes several balls located out- 
side the square acquisition surface S, which does not influence the recon- 
struction inside S (right). 



We also denote by E the operator of harmonic extension of functions from 
S to Q. I.e., for a function on S the function E<f> is harmonic inside Q and 
coincides with on S. 

Since we are dealing with the unobstructed wave propagation in the whole 
space (the surface S is not truly a boundary, but just an observation surface), 
and since we assumed that the sound speed is non-trapping and constant at 
infinity, the local energy decay type estimates of [51, 52] (see also [53, Theo- 
rem 2.104]) apply. They also lead to the following reconstruction procedures: 

Theorem 2. [14] 

1. The function f(x) in (2) can be reconstructed inside fl as follows: 

oo 

f{x) = (Eg\ t=0 ) - J A~\ sin(TA^)E(g tt ){x,T)dT. (24) 
o 

2. Function f(x) can be reconstructed inside Q from the data g in (2), as 
the following L 2 (Q)- convergent series: 

f(x) = J2fkM^), (25) 

k 
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where the Fourier coefficients f k can be recovered using one of the fol- 
lowing formulas: 

fk = K 2 9k(0) - X k 3 Jsin (\ k t)g'>(t)dt, 
o 

oo 

< fk = K 2 9k(0) + A^ 2 / cos (X k t)g' k (t)dt, or 

o 

oo oo 

fk = -K 1 f sin (\ k t)g k (t)dt = -A^ 1 / f sin (\ k t)g(x, t)^(x)dxdt, 

5 

(26) 

where 

9k(t) = J g(x,t)^-(x)dx 
s 

and n denotes the external normal to S. 

Remark 3. The function E(g a ) does not belong to the domain of the op- 
erator A. The formula (24), however, still makes sense, since the operator 
A~2 sin (tAs) is bounded in L 2 . 

This theorem in the particular case of the constant sound speed, implies 
the eigenfunction expansion procedure of [50] described in the previous sec- 
tion. However, unlike [50], it also applies to the variable speed situation and 
it does not require knowledge of a whole space Green's function. Similarly to 
the method of [50] discussed in the preceding section, this procedure yields 
correct reconstruction inside the domain, even if a part of the source lies 
outside. 

6 Partial data. "Visible" and "invisible" sin- 
gularities 

One can find a more detailed discussion of this issue for TAT in [9, 44]. Here 
we provide only a brief summary. 

Uniqueness of reconstruction does not necessarily mean the possibility of 
practical reconstruction, since the reconstruction procedure can sometimes 
be unstable. This is true, for instance, in problems of electrical impedance 
tomography, and in incomplete data problems of X-ray tomography and 
TAT [20,34,41,54]. 
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Microlocal analysis done in [10,55] (see also [56]) shows which parts of 
the wave front of a function / can be recovered from its partial X-ray or 
TAT data (see also [44] for a practical discussion). We describe this result in 
an imprecise form (see [10] for precise formulation), restricted to the case of 
jump singularities (tissue interfaces) only. 

According to [10, 55], for such singularities a part of the interface is stably 
recoverable (dubbed "visible" or "audible"), if for each point of the interface 
there exists a sphere centered at S and tangent to the interface at this point. 
Otherwise, the interface will be blurred away (even if there is a uniqueness of 
reconstruction theorem). Indeed, if all spheres of integration are transversal 
to the interface, the integration smooths the singularity, and thus reconstruc- 
tion of this interface becomes unstable. The Figure 5 shows an example of 
a reconstruction from incomplete spherical mean data. The simulated trans- 
ducers in this experiment were located along a 180° circular arc (the left 
half of a large circle surrounding the squares). In this figure the sides of 
the squares that are not touched tangentially by circles centered on S are 
noticeably blurred; any kind of de-blurring technique would not be stable in 
this context. 




Figure 5: Effect of incomplete data: the phantom (left) and its incomplete 
data reconstruction. 
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7 Range conditions 



This paper would not be complete without mentioning the intimate relation- 
ship of inversion problems with range conditions. Indeed, as it has already 
been mentioned, recovery of / from the data g is impossible, if considered as 
an inverse problem for the wave equation problem inside the cylinder S x M + . 
The possibility of inversion depends upon the fact that the solution of the 
wave equation lives in the whole space, and S is just the observation surface, 
rather than a true boundary. In other words, the data g(x,t) comes from a 
very small (infinite co-dimension) subspace in any natural function space on 
the lateral boundary S x IR + . Thus, range conditions must play a significant 
role. Indeed, they lead the authors of [14] to their results. We thus provide 
here a brief sketch of range results, following essentially the corresponding 
section of [9] . 

As it has just been mentioned, the ranges of Radon type transforms, 
including the spherical mean operator, are usually of infinite co- dimension 
in natural function spaces (in other words, ideal data should satisfy infinitely 
many consistency conditions). Information about the range is important for 
many theoretical and practical purposes (reconstruction algorithms, error 
corrections, incomplete data completion, etc.), and has attracted a lot of 
attention (e.g., [20,34,38-41,54,56-64]). 

For example, functions g from the range of the standard Radon transform 

f(x)^g(s,u) = J f(x)dx,\u\ = I, 

X-LU = S 

satisfy two types of conditions: 

1. evenness: g(—s, —uj) = g(s, uj) 

2. moment conditions: for any integer k > 0, the kth moment 

oo 

G k (u) = J s k g(u, s)ds 

— oo 

extends from the unit circle of vectors uj to a homogeneous polynomial 
of degree k in uj. 
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Although for the Radon transform the evenness condition seems to be 
"trivial", while the moment conditions seem to be the most important, this 
perception is misleading. Indeed, for more general transforms of Radon type 
it is often easier to find analogs of the moment conditions, while counterparts 
of the evenness conditions could be elusive (see [20,34,41,60,61,63]). This 
is exactly what happens with the spherical mean transform R s . 

An analog of the moment conditions was first present implicitly in [27, 
65,66] and explicitly formulated as such in [67,68]: 

Moment conditions on data g(x,r) = Rsf(x,r) in M, n are: for any 
integer k > 0, the moment 

oo 

M k (x) = J r 2k+n - 1 g(x,r)dr,xe S 
o 

can be extended from S to a (non-homogeneous) polynomial Qk(x) of degree 
at most 2k. 

These conditions are incomplete, and infinitely many others, which play 
the role of an analog of evenness, need to be added. 

Complete range description for Rs when S is a sphere in 2D was found 
in [69] and then in odd dimensions in [70] . They were then extended to any 
dimension and provided several interpretations in [26]. These conditions, 
which happen to be intimately related to PDEs and spectral theory, are 
described below. 

Let B be the unit ball in IR ra , S = dB the unit sphere, and C the cylinder 
B x [0,2] (see Fig. 6). 

Consider the spherical mean operator Rs'. 

R s f(x, t) = G(x, t) = [ f(x + ty)dA(y). 

J\v\=i 

If G(x,t) is defined by the same formula for all x G R™, then it satisfies 
Darboux (Euler-Poisson- Darboux) equation [30,71,72] 

G tt + (n- l)t- l G t = A X G. 

Inside the cylinder C, G(x, t) vanishes when t > 2 (since the spheres of 
integration do not intersect the support of the function when t > 2). 

Theorem 4. [26] The following four statements are equivalent for any func- 
tion g G C^°(S x [0, 2}), where S is a sphere: 
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Figure 6: An illustration to the range description. 



1. Function g is representable as Rsf for some f G Cq°(B). 

2. (a) The moment conditions are satisfied. 

(b) The solution G(x,t) of the interior Darboux problem satisfies the 
condition 



for any eigenfunction 4>(x) of the Dirichlet Laplacian in B. 

3. (a) The moment conditions are satisfied. 

(b) Let —A 2 be an eigenvalue of Dirichlet Laplacian in B and ip\ 
the corresponding eigenfunction. Then the following orthogonality 
condition is satisfied: 



Here j p (z) = c p -^ is the so called spherical Bess el function. 

4- (a) The moment conditions are satisfied. 

(b) Let 7)(x,X) = J g(x,t)j n /2-i(\t)t n ~ 1 dt. Then, for any m e Z, 
the m th spherical harmonic term ~g m (%, A) of g(x, A) vanishes at 
non-zero zeros of Bessel function J m+n / 2 _i(A). 



lim 




B 




(27) 



Sx[0,2] 
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One can make several important comments concerning this result (see [26] 
for a detailed discussion). In all of the remarks below, except the third one, 
the observation surface S is assumed to be a sphere. 

1. If the dimension n is odd, then conditions (b) alone suffice for the 
complete range description, and thus they imply the moment conditions 
as well. (A similar earlier result was established for a related transform 
in [70].) It is not clear at the moment whether this is holds true in even 
dimensions. 

2. The range descriptions for R s work in Sobolev scale, i.e. they describe 
the range of the operator R s : H s comp (B) ^ H!^p~ 1)/2 (S x R+). (This 
uses a recent work by Palamodov [73]). Notice that in this result it is 
assumed that the function / vanishes in a neighborhood of S, while in 
the previous theorem it was allowed for the support of / to reach all 
the way to the sphere S. 

3. If S is not a sphere, but the boundary of a bounded domain, the range 
conditions 2 and 3 of the previous Theorem are still necessary for the 
data g to belong to the range of R S - They, however, might no longer 
suffice for g to belong to the range. 

4. A different wave equation approach to the range descriptions can be 
found in [70]. 

8 Concluding remarks 
8.1 Uniqueness 

As it has already been mentioned, the uniqueness questions relevant for TAT 
applications are essentially resolved. However, the mathematical understand- 
ing of the uniqueness problem for the restricted spherical mean operators Rs 
is still unsatisfactory and open problems abound [9,27]. For instance, very 
little is known for the case of functions without compact support. The main 
known result is of [25], which describes for which values of 1 < p < oo the 
uniqueness result still holds: 

Theorem 5. [25] Let S be the boundary of a bounded domain in W 1 and 
f E L p (W l ) such that R s f = 0. If p < 2n/(n - 1), then f = (and thus S 
is injectivity set for this space). This fails for any p > 2n/(n — 1). 
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The three- and higher-dimensional uniqueness problem for non-closed ob- 
servation surface S is also still open [9,27]. 

8.2 Inversion 

Albeit closed form (backprojection type) inversion formulas are available now 
for the cases of S being a plane (and object on one side from it), cylinder, and 
a sphere, there is still some mystery surrounding this issue. For instance, it 
would be interesting to understand whether (closed form, rather than series 
expansion) backprojection type inversion formulas could be written for non- 
spherical observation surfaces S and/or in the presence of a no n- uniform 
background v s (x). The results presented in Section 1.5 seem to be the first 
step in this direction. 

The I. Gelfand's school of integral geometry has developed a powerful 
technique of the so called k operator, which provides a general approach 
to inversion and range descriptions for transforms of Radon type [39,57]. 
In particular, it has been applied to the case of integration over various 
collections ("complexes") of spheres in [39,74]. This consideration seems to 
suggest that one should not expect explicit closed form inversion formulas for 
Rs when S is a sphere. However, such formulas were discovered in [29, 45, 46]. 
This apparent controversy (still short of contradiction) has not been resolved 
completely yet. 

B. Rubin has recently discovered an alternative interesting approach to 
inversion formulas of the type of (8)- (9) for the case when S is a sphere. It 
relies upon the idea of regarding the spherical mean operator as a member 
of a broader family of operators [75] . 

In 3D, if the sound speed is constant, the Huygens' principle applies, i.e. 
the pressure p(t, x) inside S becomes equal to zero for any time T larger 
than the time required for sound to cross the domain. Thus, imposing zero 
conditions on p(t,x) and p t (t,x) at t = T and solving the wave equation 
(2) back in time with the measured data g as the boundary values, one 
recovers at t — the source f(x). This method has been implemented in 
[76]. Although in even dimensions or in presence of sound speed variations, 
Huygens' principle does not apply, one can find good approximate solutions 
by a similar approach [77]. 

A different approach to TAT inversion is suggested in [78]. It is based on 
using not only the measured data g on S x IR + , but also the normal derivative 
of the pressure p on S. Since this normal derivative is not measured, finding 
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it would require solving the exterior problem first and deriving the normal 
derivative from there. Feasibility and competitiveness of such a method for 
TAT is not clear at the moment. 

8.3 Stability 

Stability of inversion when S is a sphere surrounding the support of f(x) is 
the same as for the standard Radon transform, as the results of [9, 26, 73] 
show. However, if the support reaches outside, in spite of Theorem 1 that 
guarantees uniqueness of reconstruction, stability for some parts of f(x) lying 
outside S does not hold anymore. See [9,10,26,55] for details. 

8.4 Range 

The range conditions 2 and 3 of Theorem 4 are necessary also for non- 
spherical closed surfaces S and for functions with support outside S. They, 
however, are not expected to be sufficient, since the arising instabilities indi- 
cate that one might expect non-closed ranges in some cases. 
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